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We study the (3+l)-dimensional evolution of non-Abelian plasma instabilities in the presence of 
a longitudinally expanding background of hard particles using the discretized hard loop framework. 
The free streaming background dynamically generates a momentum-space anisotropic distribution 
which is unstable to the rapid growth of chromomagnetic and chromoelectric fields. These fields 
produce longitudinal pressure that works to isotropize the system. Extrapolating our results to 
energies probed in ultrarelativistic heavy-ion collisions we find, however, that a pressure anisotropy 
persists for a few fm/c. In addition, on time scales relevant to heavy-ion collisions we observe 
continued growth of plasma instabilities in the strongly non-Abelian regime. Finally, we find that 
the longitudinal energy spectrum is well-described by a Boltzmann distribution with increasing 
temperature at intermediate time scales. 
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I. INTRODUCTION 

One of the major outstanding questions in the theo- 
retical understanding of ultrarelativistic heavy ion colli- 
sions concerns the thermalization and isotropization of 
the quark gluon plasma. Empirical evidence in favor of 
fast thermalization and isotropization was provided by 
ideal relativistic hydrodynamical models. The success 
of these models to describe the collective flow observed 
at the Relativistic Heavy Ion Collider (RHIC) suggested 
that one generated thermal and isotropic matter at time 
scales on the order of 0.5 fm/c after the initial nuclear 
impact PHI]- Based on this success there was a con- 
certed effort to include corrections due to the finite shear 
viscosity of the plasma [5H27] . Second order viscous hy- 
drodynamics is now widely used to model collisions at 
both RHIC and the Large Hadron Collider (LHC). 

In recent years, however, studies have shown that 
there is an insensitivity to the assumed momentum space 
anisotropy of the plasma, with the data also being con- 
sistent with initially large momentum-space anisotropies 
[32 [25]. In addition, studies based on the conjectured 
anti de Sitter/conformal field theory (AdS/CFT) corre- 
spondence have shown that, although viscous hydrody- 
namical behavior emerges quickly in the strong coupling 
limit, there arc still sizable momentum-space anisotropies 
present that persist over the entire lifetime of the plasma 
[25H5T] . Based on this, extensions of viscous hydrody- 
namics that can accommodate large momentum-space 
anisotropies have been developed [2H1 !52Tf5§] . Currently 
the question of the degree of momentum-space isotropy 
of the quark gluon plasma generated in heavy ion colli- 
sions is an open question. In this paper we study the role 
played by collective unstable modes of the chromomag- 
netic and chromoelectric fields in restoring momentum- 
space isotropy of an expanding quark gluon plasma 



(QGP). 

It has been shown using both kinetic theory and dia- 
grammatic methods that when the local particle distri- 
bution function of a weakly-coupled QGP is anisotropic 
in momentum space, the system is unstable to the rapid 
growth of soft gauge fields dOHSD]- This instability has 
been dubbed the chromo-Weibel instability in reference 
to the Abelian analogue of this instability first discussed 
by Weibel [5T]. In the weak-field regime the chromo- 
Weibel instability initially causes exponential growth 
of transverse chromomagnetic and chromoelectric fields; 
however, due to non-Abelian interaction between the 
fields, exponentially growing longitudinal chromomag- 
netic and chromoelectric fields are induced that grow 
at twice the rate of the transverse field configurations. 
As a result, one finds strong gauge field self-interaction 
at late times due to high-amplitude chromoelectric and 
chromomagnetic fields and in order to reach quantitative 
conclusions numerical simulations are necessary. 

The initial numerical studies of the time evolution of 
the chromo-Weibel instability were performed assuming a 
static momentum-space anistropic (non-expanding) sys- 
tem and utilized discretizations of the gauge-invariant 
hard-loop action. The hard-loop action used includes 
the self-consistent gauge-invariant modification of all n- 
point functions in the hard- loop limit |52) . The resulting 
discretized dynamical equations were solved in tempo- 
ral axial gauge using a regular lattice to describe space 
and either a discrete lattice [55H55] or an expansion in 
spherical harmonics ">(> (i() to describe the velocity space 
of the hard particles. From the three-dimensional static 
box simulations one found that exponential field growth 
ceased when the vector potential amplitude became on 
the order of ^non-Abelian ~ Ps/ 9 ~ VfhPh, where p h is 
the characteristic momentum of the hard particles, e.g. 
Ph ~ Qs for color glass condensate (CGC) initial con- 
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ditions, fh is the angle-averaged occupancy at the hard 
scale, and p s is the characteristic soft momentum of the 
fields (p s ~ g^/JhPh)- This partial-saturation occurs at a 
scale where the chromo-fields are not yet strong enough 
to have 0(1) effects on the hard particle distribution, 
suggesting that isotropization in non-Abelian plasmas is 
parametrically slower than in the Abelian case. After 
the exponential growth ceased, a slower linear growth of 
field energy densities was observed. This linear growth 
was associated with a cascade of energy pumped into the 
soft modes to higher momentum modes through nonlin- 
ear gauge- field self-interactions [5SH55] . The resulting 
spectrum of soft gauge field excitations was shown to 
have a power law spectrum scaling like / ~ aj x p~ for 
SU(iV c ) with N c G {2,3,4,5} [Ml E3 ED- 

The presence of instabilities in weakly- coupled 
momentum-space anisotropic systems seems to be generic 
and independent of the hard-loop approximation, the 
gauge group, and, in large part, the type of theory con- 
sidered (including, of course, the weak-coupling limit of 
supersymmetric gauge theories 02 (i.'i ). They have been 
observed in numerical solutions to the full Boltzmann- 
Vlasov equations that go beyond the hard-loop approx- 
imation [B3HSS1- Analogous instabilities have been ob- 
served in numerical simulations using classical lattice so- 
lutions to pure Yang-Mills dynamics [STHTT] . A detailed 
understanding of the chromo-Weibel instability's effect 
on the isotropization and thermalization of the matter 
created in ultrarelativistic heavy ion collisions is of up- 
most importance. There have been many works that have 
addressed pieces of the puzzle (47j EH I72H75] . Recently 
there has been a highly impressive effort to parametri- 
cally estimate the effect of plasma instabilities on the 
quark gluon plasma thermalization time [7SJ [77] ; how- 
ever, being a parametric estimate it does not yet fully 
answer the question or lend itself to extrapolations to 
realistic couplings. 

In order to understand the precise role the chromo- 
Weibel instability plays in ultrarelativistic heavy ion col- 
lisions it is necessary to include the effect of the strong 
longitudinal expansion of the matter, particularly during 
its earliest stages. For the first few fm/c of the quark 
gluon plasma's lifetime the longitudinal expansion dom- 
inates the transverse expansion which only starts to be- 
come important at time scales on the order of 4-5 fm/c. 
Therefore, to good approximation, one can understand 
the early time dynamics of the quark gluon plasma by 
only considering longitudinal dynamics. The first study 
to look at the effect of longitudinal expansion was done in 
the context of pure Yang-Mills dynamics initialized with 
CGC initial conditions onto which small-amplitude ra- 
pidity fluctuations were added [781 ES] ■ The initial small- 
amplitude fluctuations result from quantum corrections 
to the classical dynamics [801 [EI]- Numerical studies have 
shown that adding spatial-rapidity fluctuations results 
in growth of chromomagnetic and chromoelectric fields 
with amplitudes ~ exp(2m° D WtJQs) where m° D is the 
initial Debye screening mass and r is the proper time. 



This growth with exp(-y/r) was predicted by Arnold et 
al. based on the fact that longitudinal expansion dilutes 
the density, thereby causing the chromo-Weibel unstable 
growth rate decrease in time [4"7] , 

Since the pioneering study of Refs. [TSJ [73] others are 
now investigating the evolution of instabilities in classical 
Yang-Mills [7IJ [5J] and scalar <\> A [53] including longitu- 
dinal expansion. In addition, a parallel effort to incorpo- 
rate longitudinal expansion into the hard-loop framework 
was begun with the first results being semi-analytic solu- 
tions for Abelian theories that also showed the character- 
istic exp(y / r) growth seen in the earlier classical Yang- 
Mills simulations, as well as rather complex early-time 
behavior [84] , In the hard- loop framework the longitu- 
dinal expansion has thus far been included only in the 
limit that the hard particles are free streaming. In this 
case it is possible to introduce a set of auxiliary variables 
similar to the static hard-loop W fields which account 
for the time-evolving momentum-space anisotropy of the 
hard particle distribution. 

The Abelian semi-analytic solutions of Ref. [84] were 
shortly followed by numerical solutions of the resulting 
coupled SU(2) Vlasov- Yang-Mills equations in the sim- 
plified case that the vector potential A and its conju- 
gate momenta II were homogeneous in the transverse di- 
rections |85| . Coupling these transversally-homogeneous 
fields to the fully three-dimensional hard-particle velocity 
distribution resulted in "1D+3V" simulations of the re- 
sulting dynamicsj^] This study found that, in the case of 
non-Abelian SU(2) fields, one also observed growth with 
exp(y / r) that was only briefly curtailed when the magni- 
tude of the transverse and longitudinal gauge field ener- 
gies became of the same order. In addition, the 1D+3V 
simulations did not see a Kolmogorov cascade at late 
times. 

The problem with such dimensionally-reduced studies 
is that they can be misleading. In fact, one finds in the 
static box case very different late time behavior if one al- 
lows for either effective one-dimensional dynamics or fully 
three-dimensional dynamics. One is therefore motivated 
to determine the full 3D+3V dynamics in the presence 
of a longitudinally expanding background. In addition, 
since the 1D+3V paper was written it was realized that 
the initial conditions used were not sufficiently generic 
and that including initial current fluctuations dramat- 
ically reduces the previously observed delayed onset of 
growth of unstable modes [HS] . One would therefore like 
to also use this type of initial condition in the full study. 

In this paper, we present the necessary 3D+3V dynam- 
ical equations for so-called hard-expanding- loops (HEL), 
discretize them in r-77-xj^ coordinates, and solve them 
numerically. For this purpose we use anisotropic lattices 
with spatial sizes on the order of N]_ x N„ ~ 40 2 x 128. 



Since, in practice, the ultrarelativistic limit |v| — > 1 is used, 
the three-dimensional velocity space is further reduced to a two- 
dimensional space (the surface of a three-sphere) . 
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At each point on the lattice we also have auxiliary fields 
W that are discretized on a velocity-lattice with size 
N<f, x N u ~ 32 x 128 amounting to 4096 auxiliary fields 
per lattice site. Needless to say this presents a com- 
putational challenge that requires parallelization of the 
resulting code. For the initial conditions we use variants 
of the initial conditions specified in Ref. [55] in which we 
have added the possibility of initializing an adjustable 
spectrum of discrete longitudinal fluctuations. As in our 
previous studies, the dimensional parameters necessary 
to fix the initial conditions such as the gluon number 
density etc. are obtained within the CGC framework. 

We find that, apart from a delay of the onset of the 
unstable mode growth due to transverse dynamics, the 
overall behavior of the three-dimensional solutions is very 
similar to the one-dimensional case. We find that the 
chromo-Weibel instability acts to restore isotropy in the 
system by inducing large longitudinal field pressure. In 
contrast to the fixed-anisotropy 3D+3V studies, we do 
not see a saturation of the instability on time scales rel- 
evant for heavy ion collisions. In order to address the 
question of the spectrum of the resulting field configura- 
tions we study the longitudinal Fourier-modes of the en- 
ergy density. We find that the longitudinal energy spec- 
trum looks like a Boltzmann distribution while remain- 
ing anisotropic in momentum space. Extrapolating to 
energies appropriate for LHC collisions, we find that the 
momentum-space anisotropy persists for approximately 
6 fm/c. We show that the isotropization time is pri- 
marily determined by the assumed magnitude of initial 
current fluctuations. 

The structure of the paper is as follows: In Section [TT] 
we briefly review the expectations one has for unstable 
mode growth in an expanding background. In Section 
|III| we review the derivation of the hard-loop equations 
of motion in a longitudinally free streaming expanding 
background. In Section |IV| we discuss the method we 
used to fix the physical scales in our simulation and dis- 
cuss the initial conditions used. In Section [V] we define 
the various observables that we will measure during the 
lattice evolution. In Section |VI| we present our main re- 
sults and interpret our findings. In Section |VII| we con- 
clude and give an outlook for the future. In three appen- 
dices we collect details concerning the numerical solution 
of the lattice equations of motion. 



II. GENERAL DISCUSSION 

Before proceeding to the presentation of the hard loop 
equations of motion and their subsequent numerical solu- 
tion, we will quickly review the presence of instabilities in 
a momentum-space anisotropic plasma and consider how 
this changes in an expanding plasma. In a longitudinal 
free streaming expansion the soft scale is time-dependent. 
Since the density of the free streaming particles drops like 
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FIG. 1. (Color online) Unstable mode growth rates (a) 
r a /m_o and (b) r_/mu for f; = 10 as a function of k z /m,D 
and 9 — arctan(fcr/fcz) where uid is the Debye mass at the 
proper time Ti so . 

n ~ 1/t and m 2 D (T) oc n/phard, we have 

f T V 1/2 

m D (T)~m D ( , (2.1) 

where mo is the "isotropic" Debye mass defined at a time 

7" ^iso ■ 

At a given proper time we can quantify the degree of 
plasma anisotropy via £ 

where px and p z are the transverse and longitudinal 



(beamline direction) particle momenta in the local refer- 
ence frame. For a longitudinal free streaming expansion 
Pt is constant while p z ~ 1/r and as a result one has 

As we will discuss in Section IV we assume that a 
plasma description becomes possible after a finite point 
in proper time To- The ratio Ti so /ro then parametrizes the 
initial momentum-space anisotropy. If this were equal to 
one, the plasma would start out isotropic and become 
anisotropic with £ > at subsequent times. However, 
motivated by the results obtained within the CGC frame- 
work [57] we consider the case that the plasma already 
has a strong oblate (£ > 0) momentum anisotropy at To, 
which will be modeled by having Tj SO <C tq regardless of 
the fact that a plasma description is certainly not possi- 
ble at times earlier than tq. By the same token, mrj, the 
isotropic Debye mass at the (fictitious) time Ti so , is just a 
parameter characterizing our free streaming background 
of hard plasma particles. 

At a given proper time r, and hence fixed plasma 
anisotropy, there is a three-dimensional band of soft un- 
stable modes associated with a fluctuation wave vector k. 
For an oblate distribution the unstable modes with the 
largest growth rate have k || n where n is the anisotropy 
direction [46]. The oblate unstable modes can be classi- 
fied as either transverse magnetic (a) or mixed (-) modes. 
The mixed modes with finite transverse momentum ex- 
tend out from the anisotropy direction to a fixed angle of 
8 = arctan(£>r/fcz) = tt/4 beyond which they are stable. 
The a-modes, on the other hand, are unstable for any 
transverse momentum. 

In Fig. [l]we show the range of unstable modes for both 
types of modes. We show the case of £ = 10 with the un- 
derstanding that the qualitative features are the same 
for all £ > 0. In a longitudinally expanding plasma lon- 
gitudinal momenta are redshifted in time, but transverse 
momenta are unaffected. As a result, the mixed unsta- 
ble modes which have any finite transverse momentum 
will eventually become stable. The a-mode growth rate 
decreases rapidly as one increases 9, so while they are 
technically unstable at all times, the growth rate of any 
mode which is not purely longitudinal becomes negligible 
at late times. Thus, at late times the system will be dom- 
inated by the dynamics of unstable modes with (nearly) 
longitudinal wave vectors]^] 

In order to gain a qualitative understanding of the dy- 
namics we can therefore focus our attention on the un- 
stable mode spectrum for purely longitudinal modes. In 



The magnitudes of px and p z stated are the "expected" values 
for the transverse and longitudinal momentum of a particle in the 

system. These can be be defined formally as px = y (p' T ) iiso/ n 

and p z = (p2) n; so /n where the averages represent integrals 
using the one-particle distribution and n is the number density. 
3 For a more detailed discussion of the dynamics of stable and 
unstable modes in an anisotropically expanding plasma see 
Ref. HI. 
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FIG. 2. (Color online) Unstable mode growth rate V /mo for 
fixed £ as a function of kz/mo where mo is the Debye mass 
at the proper time Ti so . 



Fig. [2] we plot the unstable mode growth rate for purely 
longitudinal modes for £ G {10°, 10 1 , 10 2 , 10 3 , 10 4 }. From 
this figure we can see that there is a band of modes with 
positive unstable growth rate for longitudinal momenta 
k z G (0, fc 2jmax ) and there is a well-defined maximum 
growth rate T* at each value of £. As £ increases fc 2jlnax 
increases and for £ > 10 2 one finds that T* decreases 
monotonically. This means that in an expanding plasma, 
more and more modes will become unstable as a function 
of proper time, but at the same time their growth rate is 
being reduced by the dilution of the plasma due to the 
longitudinal expansion. 

It is possible to derive asymptotic relations for k z , max 
and T* for large £ [M] . One finds that 



lim k z 



m D (l + 1/4 . 



(2.3) 



Using this we can determine the approximate proper time 
dependence of k z max for a longitudinal free streaming 
expansion 



lim k. 



- m D 

r>r iso V r iso 



1/2 



Applying the same methodology to T* one finds 

• ~ 1/2 

lim r* ~ to rj (r) nu mo — 

r»r iso \ tu 



(2.4) 



(2.5) 



As a result, we can estimate the late time unstable growth 
by integrating T* to obtain 

N(t) ~ exp (^m D jf dr' (j^J ' , 



exp (2m D ^TT iso ) , 



(2.6) 
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where we have suppressed an overall multiplicative con- 
stant. We, therefore, see that the primary effect of longi- 
tudinal expansion will be to change the late time growth 
from being a pure exponential, as was the case in a static 
box, to exp(y / T). To determine the precise nature of the 
dynamics on time scales relevant for heavy ion collisions, 
however, requires determining the full time evolution of 
all stable and unstable modes and properly taking into 
account their interactions. We will now recall the deriva- 
tion of the necessary dynamical equations to be solved 
numerically. 



III. HARD-EXPANDING-LOOP EQUATIONS 
OF MOTION 

Our study is based upon the hard-loop approximation, 
which assumes a separation of scales between the mo- 
menta of hard particles p s and the momenta of soft col- 
lective fields p s ~ gVJhPh "C Ph by a sufficiently small 
gauge coupling g. This separation obviously requires that 
fh is parametrically smaller than 1/g 2 . In an anisotropic 
plasma, fh is moreover direction dependent and what 
actually enters in the calculation of the parameters at 
the soft scale is gradients dfh/dph- In terms of the 
anisotropy parameter £ this means that at parametrically 
large £ the hard-loop approximation is applicable only as 
long as fh is parametrically smaller than 1/g 2 . 

Because we are interested in investigating within the 
hard-loop framework the earliest stages of the evolution 
of a quark-gluon plasma, which according to the CGC 
framework is born with overpopulated distribution func- 
tions and with large anisotropy, we shall treat the degree 
of anisotropy formally as being of order 1 compared to g, 
and fh of order g~ 2+e . Eventually, we boldly extrapolate 
our results to the very limits of the hard-loop framework 
by setting e = and matching with CGC parameters for 
the initial density and a strong coupling g that is nu- 
merically even larger than 1 ^] T his matching to CGC 
parameters is specified in Sec. IV in the following we re- 
capitulate the hard-expanding-loop equations, which we 
have discussed in detail before in Ref . 85J , and make the 
resulting equations explicit for the case at hand, the fully 
(3+l)-dimensional evolution. 



A. Longitudinally expanding free streaming 
background solution 

In the hard-loop approximation, the color neutral 
background distribution function / (p,x, t) for the hard 



plasma particles has to satisfy 

»-3/o(p,x,i) = 0, v»=p»/p°. (3.1) 

This is trivially solved by a stationary distribution which 
only depends on the momenta. Another solution is ob- 
tained by considering a plasma with boost-invariant lon- 
gitudinal expansion, which we take as an approximation 
for the initial stage of a heavy ion collision where the 
transverse extent of the system is taken as sufficiently 
large. Assuming isotropy in transverse directions, /o, 
which is a Lorentz scalar, can be written as [HI EH] 

/o(p,aO = fo(px,P z ,z,t) = fo(p±, P ' z ,T) (3.2) 

where the Lorentz-boosted longitudinal momentum is 

P"=7(P*-0P°), P = */t, l = t/r, r=Vt 2 -z 2 , 

(3.3) 

with p° = \Jp\ + (p z ) 2 for ultrarelativistic (massless) 
particles. 

Switching to comoving coordinates 



< = TC0sh?7, /3 = tanh?7, 
z = t sinh r), 7 = cosh 77, 



(3.4) 



we have curvilinear coordinates x a = (x T , x l , x n ) — 
(r, x , x 2 , 77) where here and elsewhere in the text indices 
i, j, . . . correspond to the two transverse spatial directions 
while Greek indices from the beginning of the alphabet 
refer to the comoving spacetime coordinates. In these 
new coordinates the metric reads 



ds 2 = dr 2 — dx 2 ! 



- T 2 dr] 2 = g a0 {T)dx a dx fi , (3.5) 



but we shall continue to write our equations explicitly in 
terms of ordinary derivatives and not deal with space- 
time covariant derivatives. The gauge covariant deriva- 
tive thus always mean^]z? Q = d a — ig[A a , ■]. 

The field strength tensor is defined as F a p = d a Ap — 
dpA a — ig[A a , Ap] also in the comoving coordinates (with 
all indices down), in which the non-Abelian Maxwell 
equations can be written compactly as 



D a (rF^)=f, 



(3.6) 



where indices have been raised with the inverse of the 



metric g a p(r) introduced in Eq. (3.5|. 

Similarly to space-time rapidity 77, we define momen- 
tum space rapidity y for the massless particles according 
to 



p M = p± (cosh y, cos 4>, sin <j>, sinh y) . 



(3.7) 



4 In the notation of Ref. [77] where / ~ a~ c , r <x aj a , ~ 
= aj d , our framework is located at parametric time a = 
with parametric occupancy c = 1 — | and parametric anisotropy 
d = 0. 



5 The relation to 3- vectors is defined by d a = d/dx a and = 
(<f>,Z). Thus A a = (A T ,-A x ,—A v ,An). 
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In comoving coordinates, this reads 



if it satisfies 



— coshryp — sinh r]p z = p± cosh(y — rf), (3. 
p 11 = — Pti/t 2 = (cosh r\ p z — sinh?7p )/r 



p' z /t — p± sinh(y — rf)/n 



(3.9) 



Instead of the standard light-like vector = p^/p 
which contains a unit 3-vector and which was introduced 



in Eq. (3.1), we shall define 



P a ( 1 
V a = — = cosh(y — rf), cos <^>, sin</>, — sinh(?/ — 77) 

p± V T 

(3.10) 

normalized such that it has a unit 2-vector in the trans- 
verse plane. 
Since 



P T d T Pn{x) 



= —p\ sinh(y — rf) cosh(y — rf) , 

= -pid vPr ,(x) , (3.11) 



this can be solved by / Q (p,x,£) = f {p±,p v (x)) = 
fo(p±, —p' z (x)t(x)). For the case of longitudinal free 
streaming which is isotropic at the particular proper time 
t = T iso one can write /o in the form 



fo(p,x) = fa 



ipl + C—y 



hsoi^Jpi+plh 



(3.12) 

Note that fo above falls into the general Romatschkc- 
Strickland form for momentum-space anisotropic distri- 
bution functions 1351. 



B. Gauge-covariant Boltzmann-Vlasov equations in 
a longitudinally expanding plasma 

In comoving coordinates the gauge-covariant 
Boltzmann-Vlasov equations for colored perturbations 
Sf a of a neutral collisionless plasma with boost-invariant 
background distribution fo read 



V-DSf a \ 



(3.13) 



Here the derivative on the left-hand side has to be 
taken at fixed Cartesian p M rather than fixed comov- 
ing p a . Notice also that only derivatives of fo(p±,p v ) 
with drp\ where the 4-index is up do not introduce ex- 



plicit r dependence so that one still has p ■ d (<9^/o)| p 

P -d(dU )\p = o- 



Eq. (3.13) can be solved in terms of an auxiliary field 
Wp(x\ (j>, y) that does not depend on the hard scale p° 
and which is defined by 

5f(x;p) = -gWf3{x;<j),y)d^ p) fo{pj_,p,f), (3.14) 



(3.15) 



Since the fluctuations 5f a give the induced current in 



fa=9tR 



d 3 p p 



(2tt) 3 2p< 



Jf a (p,x,t), (3.16) 



j can be expressed in terms of integrals over the W fields. 
(Here tn is a suitably normalized group factor, while the 
total number of degrees of freedom of the hard parti- 
cles is contained in the normalization of the distribution 
function /p.) 

With (pU2j) we have 



r 2 

ISO 



9f P )/o = f^pl+pyri 

I 0, — cos <p, — sin <fi, — 5~ sinh(y — rf) 
= fo S 



1 + -£r sinh 2 (y - 77) 



(3.17) 



which yields 

r = - 



m D 


r-2-K 
/{) 


d<b r°° 

^ Loo 


dyV a 


(- 


r 2 

T? 

ISO 


sinh 2 (y 





where 



W = VWi 2~^) W v , 

^"iso 

V 1 = (cos <f>, sin </>) , V v = —t sinh(y — rj) , (3.19) 



and 



ral = -g 2 t R J™^fU P )- (3-20) 

The (constant) mass parameter mo equals the Debye 
mass at the proper time Ti so . 

The combination W introduced above satisfies 



V ■ DW = \ V" l Fi 

+V l V r >F lv [ 1 



„ -r -^V^F nT ) V T 



(3.21) 



This single equation for W in combination with the Yang- 
Mills equations and the integral giving j in terms of W 
closes our equations of motion. To solve them numer- 
ically, we adopt the comoving temporal gauge A T = 
and introduce canonical conjugate field momenta for the 
remaining gauge fields according to 



and 



W = rd T A t = -rdrA 1 = -H i; , (3.22) 

(3.23) 



IP = -d T A„ . 

T 
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In terms of fields and conjugate momenta, the Yang- 
Mills equations take the form 

rd T W = j n - D,F\ , 



(3.24) 

r^firlli = f - DjF ji - D v F vi , (3.25) 
while the Gauss law constraint takes the form 



T3 



D V BP + Dm, 



(3.26) 



In temporal gauge, where Fj T = ILj/r and F VT 
-rW, the field equation for W, Eq. (3.211, becomes 



<9 t W(t, xx, m ;<f>,y) 



1 



cosh(y) 

T 

tanh(y) 



v l DiW 



sinh(y) 



A,W 



t 2 sinh(y) 



IF 



, ,v l F lv , (3.27) 



with y = y — rj. 



In the limit that all fields are independent of the trans- 
verse spatial directions Eqs. ( 3.24 )-( 3.27) reduce to the 



1D+3V equations of Ref. [85]. 

We can recast (3.27) into a form which is more conve- 



nient for computing the currents in Eq. (3.18) by defining 



W(r, x ± , ??; 0, y) = /(r, r iso , y) W(r, x_l, V, <i>, v) , (3-28) 
with 

r 2 



f(r, r iso , y) = ( 1 + 32" sinh 2 y 



(3.29) 



We also replaced ybyy = y—r\ as argument of W because 
the auxiliary fields turn out to be peaked around y ~ t). 
Now using 



OT 



(3.30) 



D v W{t, xx, ??; 0, y) = (D, - dy) [fW(r, xx, r?; y)] 

9/ 



/ (A, - d v )W - ^_ W (3.31) 



AW = /AW , 



(3.32) 



together with 



tanh y 



df = dl 
dy dr 



we obtain 

5 r W(r,xx,»7;0,y) = - 



1 

coshy 
sinhy 

T 



t/AW 



+ 



f{r,n ao ,y) 



+ 



I^n, - r2si 2 D ^ n" 

T T*. 

ISO 

tanhy 



v'Fi, 



(3.33) 



In terms of W(r,]t±,r);<p,y) the expression for the cur- 
rent ( 3.18[ ) simplifies to 

j Q (r,xx,r/) = 

/ dyl/ Q W(r,xx,?7;^y). (3.34) 

The equations of motion listed above are numerically 
solved by discretizing them in space and velocity space 
(hence the designation 3D+3V). The gauge fields live on 
the 3-dimensional space parametrized by space-time ra- 
pidity rj and two transverse coordinates xx- The W 
field lives additionally in velocity space, which because 
of the masslessness of the hard particles is, in the end, 
2-dimensional, parametrized by y and 4>. 

For the details of the lattice discretizations used we 
refer the reader to Appendix [A] 



IV. INITIAL CONDITIONS 

A. Matching of the Debye mass with CGC 
parameters 

As in our 1D+3V simulations [85], we evolve from an 
initial time r ~ Qj 1 an d fix the density of our initial 
plasma such that it matches estimates obtained from the 
CGC framework. 

According to Ref. [90], the initial hard-gluon density 
can be written as 



"(to) 



NgQl 



(4.1) 



with c being the gluon liberation factor, which following 
an analytical estimate by Kovchegov [5T] we choose as 
c = 2 In 2 R3 1.386. While being significantly higher than 
the original estimates c ~ 0.5 of Ref. [9"2l 193] . this value 
is in fact rather close to the most recent numerical result 
c ~ 1.1 by Lappi [M] , 

In our effective field equations, the initial hard-gluon 
density enters only through the mass parameter tod, 
which is defined as the Debye mass at the proper time 
r iso . In the glasma phase of the CGC framework, the 
pressure at early times is strongly anisotropic, with the 
longitudinal pressure starting out even with negative val- 
ues. To model this approximately, we formally choose 
T iso tq, so that our initial particle distribution has ini- 
tial pressure Pj, <C Pt- Sticking to our previous choice 
in Ref. |85j we take T{ so = 0.1 To- The correspondingly 
oblate distribution function is taken to be obtained from 
hofr) = N{2N g )/{e^ T - 1), where N g = N 2 - 1 is 
the number of gluons, since in CGC calculations an ap- 
proximately thermal distribution was obtained for the 
gluon distribution in transverse directions. Following 
Ref. [55] we set this t rans verse temperature T = Q s /d 
with dr 1 ~ 0.47. Eq. (4.1) then fixes the normalization 
factor M through 



n(T ) 



To 



n(r is 



IT 2 



(4.2) 
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In a plasma containing only gluons with distribution 
function /i so , the Debye mass is given by 



m 2 D (T iso ) = N 



4?ra s N c T 2 



(4.3) 



With N c = 3 and the above values for c and d we thus 
obtain 



m 2 D (Ti so )TQ(Q s T Q y 



ncd 



6C(3) r iE 



1.285 



(4.4) 



In our previous studies of a stationary anisotropic plasma 
we have observed little difference between simulations us- 
ing gauge group SU(2) versus SU(3) provided the s ame 
value of m_o was used [SUES], so we adopt the value (4.4) 
for our simulations with gauge group SU(2). 

Notice that in the above matching which involved an 
overpopulated distribution function n(ro) oc aj 1 the 
gauge coupling dropped out in the mass parameter mf^. 
As discussed in Section |III[ this means that we are ex- 
trapolating the hard loop framework, which assumes a 
parametric separation of hard and soft scales, to its very 
limits. In the following we shall compare hard and soft 
contributions to the pressure and find that the soft field 
contributions are small compared to the hard particle 
contributions even after plasma instabilities have grown 
nonperturbatively strong. As long as this is the case, we 
assume that the hard loop framework is still applicable. 

In order to compare soft and hard contributions, we 
finally have to fix the gauge coupling. For that purpose 
we shall choose a s = 0.3 or g — 1.94 as a representative 
value. 



B. Initial field fluctuations 

In order to have seed fields for the unstable modes in 
an anisotropic plasma with oblate anisotropy, initial fluc- 
tuations that break perfect boost invariance are required. 
Fluctuations in the sources of heavy-ion collisions as well 
as vacuum fluctuations in all fields are inevitable, and 
by "natural selection" those fluctuations which lead to 
the most rapid onset of growth will dominate all later 
dynamics. 

In previous hard-loop lattice simulations with fixed 
anisotropy the question of which initial conditions to 
choose was rather unimportant as long as unstable modes 
were excited. Seed fields in chromo-fields or in W fields 
were considered on the basis of convenience. 

As it turns out, more care is needed in the expand- 
ing case. In [84] . where the formalism of hard expand- 
ing loops was introduced and studied semi-analytically 
in the (1+1) -dimensional Abelian case, only initial con- 
ditions formulated in terms of transverse electric fields 
were considered. Likewise, only seed fields in transverse 
chromo-fields were subsequently employed in the numer- 
ical 1D+3V non- Abelian lattice study of Ref. [85 , which 
in the weak-field regime reproduced the earlier semi- 
analytical results, and thus also the original finding of a 



(with regard to heavy-ion collisions) uncomfortably long 
delay of the onset of growth of plasma instabilities. (The 
generalization considered in Ref. [55], namely to also in- 
tialize magnetic fields did not change this conclusion.) 

In Ref. [86] the semi-analytical treatment of Ref. [84] 
was generalized to the much more complex case of 
generic (3+l)-dimensional Abelian modes in an expand- 
ing plasma, and at this occasion also the most general 
initial conditions were considered, involving both elec- 
tric and magnetic fields as well as the auxiliary W fields 
which describe fluctuations in the induced currents. Sur- 
prisingly enough, initial fluctuations in the W fields lead 
to a drastic (order-of-magnitude) reduction of the initial 
delay of the onset of growth. Evidently, initial condi- 
tions in the electric and magnetic fields predominantly 
give stable plasmon modes and less strongly excite the 
unstable modes. The latter are instead more easily trig- 
gered by fluctuations in the induced currents described 
by the W fields. 

The simplest initial conditions that provide seed fields 
for Weibel instabilities while having initial vanishing 
charge density are <fi and y-independent fluctuations of 
the component fields Wi(x; (p, y) and W 71 (x;(p,y). The 
former induce transverse currents which are most directly 
related to the a modes, whereas a </> and y-independent 
W v seeds longitudinal currents that are less important 
for the plasma instabilities. Because of their subdomi- 
nant effect, we have mostly omitted W„ seeds and only 
kept Wi(x;4>,y) when assembling the initial W field. 

Another point to consider is the spectrum of initial 
fluctuations. Because we are using highly anisotropic 
lattices with particularly fine resolution in the longitu- 
dinal direction, initializing with white noise fluctuations 
would correspond to very high UV noise in longitudinal 
wave numbers. We have therefore implemented an ad- 
justable mode number cutoff, A„, in wave numbers v dual 
to the rapidity variable r\ and populate all modes oc e lvri 
equally below this cutoff, with white noise in transverse 
directions. Because the "natural selection" of plasma in- 
stabilities quickly picks out the most strongly growing 
modes, we have refrained from attempts to model the 
initial spectrum other than ensuring that a good range 
of seeds is available. 



OBSERVABLES 



Here we list the quantities which we will present in the 
results section. We present only the continuum formu- 
lae. For the details of the lattice discretizations used we 
refer the reader to Appendix [A] Note that in most of the 
results presented we have averaged observables over a set 
of runs in order to account for variations in the random 
initial conditions employed. 
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A. Field energy densities and pressures 



C. Energy spectra 



The transverse/longitudinal electric and magnetic 
components of the field energy density are given by 



£ — £t + £l — £b 
= tr 



Si 



t~ 2 F 2 . 

7]l 



£r 



F 



XIJ 



(5.1) 



and the Hamiltonian density is given by H = t£. The 
transverse and longitudinal field pressures are obtained 
via 



Vl cld = £ t -£l, 

afield 



TM1C1Q c 



(5.2) 
(5.3) 



Note that from the above one has at all times 2'Pf. cld 



;>ficld _ 

traceless 



•pheid _ g such that the energy momentum tensor is 



B. Particle Pressures 

In a comoving frame, the energy density and pressure 
components of the hard particle background can be de- 
termined by evaluating 



T. 



Ct0 



part. 



(27T) 



d< PT dyp a pPf , 



(5.4) 



which yields 



£ part, (t) = T 



TT 

part. 



1 arcsin y/l — r 2 
^2 ~! 



1 



f. 



1. 



p part. (T)= _ T , art 



4(f 2 



arcsin 



2(f 2 - 1) 



1 arcsin vl — r 2 



1 



(5.5) 

f. 

(5.6) 
f iso , (5.7) 



where £i SO = £ pa rt. ("Hso), r = t /t- iso and we have assumed 
f > 1. 

In the results section as a measure of isotropization we 
will present plots of the ratio 



-jficld 



■jpart. 



(5.8) 



If this quantity is less that one, then the system possesses 
an overall oblate momentum-space anisotropy and if it is 
greater than one, then it possesses a prolate momentum- 
space anisotropy. 



In order extract spectral information about the field 
configurations, the canonical way to proceed is to gauge 
fix to a spatially smooth gauge such as Coulomb gauge 
and then extract mode occupation numbers from either 
the electric or magnetic fields 58, 96, 97J. However, such 
a method is not free from ambiguity in the infrared due 
to the lingering problem of large gauge transformations 
(aka Gribov copies). 

Here we follow a different method introduced by 
Fukushima and Gelis [52] in which we extract the elec- 
tric and magnetic fields at a given proper time from the 
lattice simulation using 



Ei{* T ,n) = T~ 1 u i , 
e l ( Xt , v ) = w, 

B x (x. t ,jj) = F n , 



iga n i 



■tr[t a (l-U vy )} 



Sy(x T ,I)) = F qx 



ti[t a (l - U vx )] , 

B L (x T , V ) = F xy ~ -tr[t Q (l - U xy )] . 
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(5.9) 



We then perform a three-dimensional Fourier transform 
of each field component, e.g. 

Ei(k T ,v)= J l^g^xx^K 1 ™ e 4 "". (5.10) 

Since we are primarily interested in the longitudinal spec- 
tra, we integrate over the transverse wavevectors to ob- 
tain, e.g. 



Ei{v) 



Having obtained the field components we can decom- 
pose the energy density in terms of the longitudinal 
wavenumber 



£e = 
£b = 



dv 
2^ 
dv 



£e(v) 



dv 
2^ 

2ir V ' / 2tt 



[Se l (v) +S Et (v)] . 
[SbM + SbM] 



(5.11) 



where we have the energy density at each longitudinal 
wavenumber 

£ El (v) = tx[E L {-v)E L {v)\ = tr\E L \ 2 , 

S Et {v)= tr[E i (-v)E i (v)]= tr |£T> 

i£{x,y} i£{x,y} 

£ Bl {v) = tv[B L (-v)B L (v)} = tr\B L \ 2 , 

S Bt (v)= ]T te[B l (-u)B i (y)]= J] tr|i?f, 



iE{x,y} 



(5.12) 
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where the traces are color traces. The total longitudinal 
energy spectra are obtained by summing all components 

£(u) = £ El (j/) + £ Et (y) + £ Bl {y) + £ Bt (v) . (5.13) 



The spectral decomposition (5.12) is not gauge invari- 



ant; gauge transformations could in principle still redis- 
tribute the energy distribution in v, but this redistri- 
bution is limited by the fact that the integrals (5.11) is 



gauge invariant. We thus expect that the degree of gauge 
dependence is much milder than in bare mode occupation 
numbers of the gauge fields before they are made maxi- 
mally smooth by going to Coulomb gauge. 

Note that one can compute the total energy density via 
Eq. (5.1) and Eq. (5.11) and compare as a crosscheck of 
the spectra calculation. Numerically we find very good 
agreement between the two methods. We have also per- 
formed a Fourier analysis of the spatial distribution of 
the (gauge-invariant) chromo-field energy on the lattice. 
Besides the expected peak at zero momentum, we found 
that the remaining spatial fluctuations reflect closely the 
spectral decomposition defined through Eqs. (5.12). 



VI. RESULTS 

In this section we present the results of our numer- 
ical simulations for SU(2) gauge fields which include: 
real-time gauge field energy densities, particle and field 
pressures, energy spectra, and fit to the energy spec- 
tra. For all results shown in this section we initialize 
current fluctuations (via VV fields) with an amplitude A 
as described in Section [IV] and Appendix [C] In order 
to generate occupation numbers ~ 1/2 consistent with 
those expected from initial quantum-mechanical rapidity 
fluctuations [5U] one should choose A ~ 1.6. Unfortu- 
nately, due to numerical limitations stemming from the 
fact that we simulate compact gauge groups, we are un- 
able to use such a large value of A. Instead in the main 
plots shown below we use an initial current fluctuation 
amplitude of A = 0.8 which can be expected to result 
in longer isotropization times than one would would ob- 
tain with the larger seed values necessary. In order to 
assess the dependence of our results on A we present the 
variation of the energy density and pressure ratio. In the 
conclusions we will discuss the extrapolation of our result 
to realistic values of A. 

For all results shown the lattice spatial size was iVf. x 
N v = 40 2 x 128 with transverse lattice spacing of a = 
Qj 1 and longitudinal lattice spacing of a„ — 0.025. The 
lattice size in velocity space was N u x N$ = 128 x 32. The 
longitudinal spectral cutoff for the current-based rapidity 



fluctuations was taken to be A v v n 



15.7. 



The initial time was taken to be r = Q~ x and we used 
Ti S o/To = 0.1. For the temporal time step we use e = 
10~ 2 To. For details of the lattice discretizations used for 
the equations of motion we refer the reader to Appendix 
|A") When plotting observables as a function of time we 
will plot them as a function of f = Q s r/10. For LHC 




) 1 2 3 4 5 6 7 
X = Q s X / 10 

FIG. 3. (Color online) Chr omoelectric, chromomagnetic, and 
total energy densities (5.1 1 as a function of proper time from 



averaging over our standard set of runs. Proper time is nor- 
malized such that when using Q„ = 2 GeV each unit of Af is 
1 fm/c. See text for simulation parameters used. 



one has Q s ~ 2 GeV = (0.1 fm) 1 and for RHIC one has 
Q s ~ 1.4 GeV = (0.14 fm)- 1 . The division by a factor of 
10 makes it so that when considering LHC energies each 
interval of Af = 1 is 1 fm/c. At RHIC each interval of 
Af = 1 is 1.4 fm/c. 

For numerical tests such as varying the lattice spacing, 
lattice size, spectral cutoffs, and velocity resolution we 
refer the reader to App. [D] The lattice equations of mo- 
tion are written in terms of rescaled dimensionless fields. 
When comparing pressures in soft fields with pressures 
from hard particles, we have assumed a value of g = 1.94 
consistent with a s — 0.3 which is in the right ball park for 
RHIC and LHC heavy ion collisions. Note that formally 
our results are only trustable in the weak-coupling limit 
and we are making a bold extrapolation when we assume 
a s = 0.3. Nevertheless, we do this in order to obtain a 
rough estimate of the isotropization time associated with 
the chromo-Weibel instability in a background which is 
undergoing longitudinal free-streaming expansion. 



A. Energy densities 

In Fig. [3] we show the chromoelectric, chromomagnetic, 
and total energy densities (5.1) as a function of proper 



time. The results shown are averaged over 50 runs which 
will serve as our standard set of runs for most observ- 
ables in this sectionj^] From Fig. [3] we see that for the 



In App. [d] Fig. |l2(a)| we plot the total field energy density re- 
sulting from all 50 runs for comparison. 
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FIG. 4. (Color online) Total field energy density for different 
initial current fluctuation magnitudes A S {0.1,0.2,0.4,0.8}. 
See text for simulation parameters used. 



FIG. 5. (Color online) Hard particle and field pressures scaled 
by Tp t as a function of proper time. The data were taken from 
the same set of runs as Fig. [3] 



first f < 1.2 the soft fields are depleted by the longi- 
tudinal expansion. After this time the unstable modes 
present in the initial condition begin to show appreciable 
growth. Initially all components of the chromofield start 
out with approximately equal energy density, but at this 
time the system begins to be dominated by transverse 
chromomagnetic fields. However, due to the large ampli- 
tude of the initial current fluctuations we quickly see the 
development of large transverse chromelectric fields fol- 
lowed by rapid growth in the longitudinal chromoelectric 
and chromomagnetic fields. 

All field components become approximately the same 
magnitude at a time of f ~ 3.5 when A = 0.8. We will 
refer to the point in time at which all components of the 
field energy density give approximately the same contri- 
bution as the "non-Abelian point" . From this point on, in 
contrast to the fixed-anisotropy simulations, one does not 
see a saturation of the exponential growth, just a mod- 
erate reduction of the growth rate. Instead we see that, 
similar to the 1D+3V simulations, the transverse chro- 
moelectric and chromomagnetic fields begin to dominate 
the energy density and do so for the rest of the simula- 
tion. As we will see below, by the end of the simulation 
a large portion of the energy is in ultraviolet longitudi- 
nal lattice modes and one starts to see lattice artifacts; 
however, up to this point we see no sign of saturation of 
the roughly exponential growth in the chromofields. 

In Fig. [4] we show the total field energy density 
for different initial current fluctuation amplitudes A € 
{0.1,0.2,0.4,0.8}. As can be seen from this figure, apart 
from a slight reduction in unstable mode growth when 
the fields reach the non-Abelian point (which moves to 
large times for smaller A), the behavior is qualitatively 
independent of the assumed amplitude. We note that 
there is a fundamental limit on how large one can make 



A without violating the assumptions of the hard-loop ef- 
fective theory we employ. In practice, this limit is set by 
the physical requirement that the majority of the energy 
density should still be contained in the hard particle dis- 
tribution function. We note that for A = 0.8 we are still 
safely below this bound with the initially induced soft 
fields only carrying only ~ 1.5% of the total energy with 
the vast majority of the energy coming from the hard 
sector. 



B. Pressures 

In Fig. [5] we show the hard particle and field pressures 
scaled by tjt as a function of proper time. The data were 
taken from the same set of runs as Fig. [3] and the pres- 
sures were computed using Eq. (5.3). The scaling chosen 



in this figure renders the vertical axis dimensionless and 
has the added benefit of making the scaled hard particle 
transverse pressure constant for better visualization. 

As can be seen from Fig.[5]the system is initially highly 
anisotropic with the transverse particle pressure domi- 
nating all other contributions. The r-scaled longitudi- 
nal particle pressure drops like 1/r 2 . Note that at early 
times the field component of the longitudinal pressure 
can become negative as evidenced by Fig. [5j This is con- 
sistent with the finding of others [52] and is a result of 
coherent field modes. Without the unstable field growth, 
the system would continue to become more and more 
anisotropic as time progresses and continue to experience 
positive and negative pressure oscillations. However, as 
Fig. [5] demonstrates, unstable field modes begin to gen- 
erate a growing longitudinal field pressure that at late 
times dominates all other pressure components. 

It should be noted, however, that by the time the Ion- 
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FIG. 6. (Color online) Total longitudinal pressure over the 
total transverse pressure as a function of proper time. The 
data were taken from the same set of runs as Fig. [3] 



gitudinal field pressure becomes of the same magnitude 
as the transverse particle pressure one already expects 
to see a significant amount of backreaction of the hard 
particles on the unstable chromofields. Physically this 
should result in a saturation of the field pressure growth 
due to energy conservation. In addition, the back reac- 
tion would serve to isotropize the particle sector. Such a 
physical saturation is, unfortunately, not describable in 
the hard-loop framework since in this framework the hard 
particles act as an energy reservoir that can continue to 
pump energy into the soft sector indefinitely. Sans this 
caveat, we believe that this result shows evidence that 
the chromo-Weibel instability can restore isotropy on the 
fm/c time scale. 

In Fig. [6] we show the total longitudinal pressure over 
the total transverse pressure (5.8) as a function of proper 
time. The data were taken from the same set of runs as 
Fig. [3] This plot condenses the information seen in the 



previous plot allowing one to easily see the point at which 
the plasma becomes isotropic in momentum space. As 
can be seen from this figure this occurs at approximately 
f = 6.5; however, the system continues to evolve beyond 
this point with the total longitudinal pressure then ex- 
ceeding the total transverse pressure. This is most def- 
initely an artifact due to the lack of the back reaction 
of the hard particles on the chromofields. Therefore, we 
are only fully confident in the results we obtain at earlier 
times. 

In Fig. [7] we show the tota l lon gitudinal pressure over 
the total transverse pressure (5.8) as a function of proper 



time for different initial current fluctuation magnitudes 
A £ {0.1,0.2,0.4,0.8}. The data were taken from the 
same runs as shown in Fig. |4j The purpose of this figure 
is to show that the variable which has the biggest effect 
on the isotropization time is our assumed magnitude of 



FIG. 7. (Color online) Total longitudinal pressure over 
the total transverse pressure as a function of proper time 
for different initial current fluctuation magnitudes A € 
{0.1,0.2,0.4,0.8}. The data were taken from the same runs 
as shown in Fig. [4] 



the initial current fluctuations, A. From this figure we 
see that each time A is doubled the time at which the 
system becomes isotropic decreases by f ~ 1. In the 
range of A's studied, the isotropization time varies loga- 
rithmically like a — b log A with a and b being constants. 



C. Energy spectra 

In Fig. § we show the run-averaged longitudinal en- 
ergy spectra obtained via (5.13) at different proper times 



as a function of (a) the longitudinal wavenumber v and 
(b) the longitudinal momentum k z = v/t. The data for 
both plots were taken from the same set of runs as Fig. [3] 
In both figures the vertical axis is logarithmic while the 
horizontal axis is linear. From Fig. 8] (a) we see the rapid 
emergence of an exponential distribution of longitudinal 
energy. The exponential spectra persist during the entire 
evolution. In Fig. [8] (b) we show the spectra as function 
of the physical momentum so that one can now see the 
effect of the red-shifting of the longitudinal momentum 
with time. In addition, from this figure we can easily 
determine a kind of effective longitudinal temperature 
which can be extracted from the slopes of the curves. 
Below we will define a fit function and extract the longi- 
tudinal temperature as a function of proper time. 

Note that the emergence of this exponential spectrum 
is not solely due to the widening unstable mode band. 
Instead having nonlinear mode-mode coupling is vitally 
important in order to populate high momentum modes 
which are rapidly becoming unstable as time progresses. 
In order to illustrate this point in Fig. [9] we show the 
corresponding spectra from Abelian runs. The lattice size 
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FIG. 8. (Color online) The longitudinal energy spectra at 
various proper times as a function of (a) v and (b) k z = v/t. 
Data taken from the averaged runs shown in Fig. [3] 



for these Abelian runs were exactly the same as for the 
corresponding non- Abelian run shown in Fig. [8} however, 
we chose a smaller value of A in order to eliminate the 
possibility of artificial nonlinearities due to the fact that 
we are simulating compact U(l). As we can see from 
this figure, only modes present in the initial conditions 
are amplified in the Abelian case, hence demonstrating 
that the emergence of an exponential longitudinal energy 
spectrum is intrinsically non- Abelian (nonlinear). 

At first sight our exponential distribution of longitu- 
dinal energy seems to be different than the result ob- 
tained by Fukushima and Gelis who saw the emergence 
of a power-law spectrum in Yang-Mills solutions in an ex- 
panding QGP 82 ; however, we note, importantly, that 
they saw the emergence of a power-law longitudinal en- 
ergy spectrum only at extremely late times correspond- 
ing to t > 150. At early times, their spectra also appear 
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FIG. 9. (Color online) The longitudinal energy spectra at 
various proper times as a function of v for Abelian runs. For 
this figure the spectra from 40 runs were averaged. 



consistent with an exponential distribution of longitudi- 
nal energy. Since we do not include the back reaction, 
we are unable to comment on the asymptotic behavior 
of the spectra since we currently see no evidence of soft- 
scale saturation of the unstable mode growth. In addi- 
tion, power law scaling usually emerges in the infrared 
and, in that sense, we are limited due to small lattices. 

In Fig. 10 we show fits to spectra shown in Fig. [8] (b) 
at several different proper times. For the fit function 
we assumed that the spectra corresponded to the energy 
density obtained from a massless Boltzmann distribution 
that has been integrated over transverse momenta 



foc^ dk z d 2 k T \jk\ + k 2 z cxp (-y/k% + k*/Tj , 

oc J dk z (k 2 z +2\k z \T + 2T 2 )exp(-\k z \/T) . (6.1) 

The integrand in the above expression was taken as our 
fit function 

£fit(*!») = A (k 2 z + 2\k z \T + 2T 2 ) exp (-\k z \/T) , (6.2) 

where we have allowed for an overall multiplicative con- 
stant A. At each proper time we fit the two parameters 
A and T; however, at early times we manually exclude 
regions of the spectra that are part of the "noisy plateau" 
at high longitudinal momenta, e.g. k z > 8 Q s from the 
f = 0.3 panel shown in Fig. [10] are excluded from the fit 
data. 

As can be seen from Fig. 10 we see evidence of a 
very rapid emergence of a Boltzmann longitudinal en- 
ergy spectrum. At f = 0.3 the fit is already working 
quite well with the bumps seen in the spectra being non- 
linear resonance "copies" of the initial theta-function-likc 
distribution of longitudinal energy. By f = 2.3 virtually 
all information about the initial condition is gone and by 
f = 3.3 the system seems to exhibit an exceptional degree 
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FIG. 11. (Color online) The time dependent longitudinal tem- 
perature extracted from the data contained in Fig.[8]using the 
fit function |6^2}. 



gins to grow, the soft sector begins to heat up. We 
note in this context that the hard particle distribution 
is highly anisotropic, making it hard to associate a tem- 
perature with. The transverse temperature given by 
Phard ~ Qs is a constant for longitudinal free streaming; 
however, one can associate a kind of isotropic temper- 
ature by computing the fourth root of the energy den- 
sity £ = 7£(£)£iso(Phard) [31] ■ One finds at late times 
(r > T iso ) that £ ~ r" 1 so that T eff , hard - £ 1 / 4 ~ t" 1 / 4 
which decreases less quickly than ideal hydrodynami- 
cal behavior for which one has T ~ r -1 / 3 . Since the 
hard particles still dominate the energy density, the com- 
bined soft plus hard effective temperature still decreases 
in time. 



VII. CONCLUSIONS 



of longitudinal thcrmalization with all information about 
the initial condition lost. We only show six specific times 
in Fig. 10 however, at all simulation times f > 0.3 the 



fits seem to work remarkably well. We note, importantly, 
that although the spectra shown in Fig. [lO] are averaged 
over runs, one sees the emergence of such a Boltzmann 
spectrum on a run-by-run basis. We have averaged over 
runs in order to remove statistical noise and improve the 
quality of the fits. 



In Fig. 11 we show the extracted fit temperatures us- 
ing (6.2) as a function of proper-time. We see from this 



figure that at early times the soft sector cools down due 
to longitudinal expansion, but once the instability be- 



in this paper we have studied the dynamics of the 
chromo-Weibel plasma instability in a longitudinally ex- 
panding plasma by numerically solving the full 3D+3V 
realtime evolution of the hard-loop equations of motion. 
We utilized current fluctuations as the initial condition so 
that the initial fields were self-consistently triggered by 
the hard particles. We had three important findings: (1) 
there is no saturation of the chromo-Weibel instability 
at the "soft-scale" on timescales relevant for heavy-ion 
collisions, (2) the dominant transverse chromomagnetic 
fields generate a rapidly growing longitudinal pressure 
that works to isotropize the system on timescales relevant 
for heavy-ion collisions, and (3) in the process of evolu- 
tion the longitudinal energy spectrum shows no signs of a 
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power-law spectrum associated with Kolmogorov turbu- 
lence, but instead shows evidence for rapid longitudinal 
thermalization of the gauge fields. 

The finding that there is no soft-scale saturation of 
the plasma instability is important since this means that 
on the time scales relevant for heavy ion collisions the 
back reaction of the hard degrees of freedom could be 
important. This suggests that it will be of the upmost 
importance to make an in depth study of the dynamics of 
an unstable expanding plasma using classical Yang-Mills 
and Boltzmann-Vlasov simulations. However, care will 
have to be taken to make sure that these simulations can 
properly describe the soft collective modes of the system 
consistent with hard-loop dynamics in the high temper- 
ature limit. The fact that we do not witness soft-scale 
saturation of the chromo-Weibel instability is consistent 
with previous analyses of plasmas possessing a fixed high- 
magnitude momentum-space anisotropy 59, 60 . In the 
case of HEL, the red shifting of the longitudinal momen- 
tum causes transverse unstable modes to become more 
and more stable as a function of time, while purely lon- 
gitudinal modes continue to grow. We cannot rule out a 
very late time saturation on timescales far beyond what 
we have studied; however, such large time scales are prob- 
ably not relevant to understanding thermalization of a 
QGP generated in heavy ion collisions. 

Our second finding concerned plasma isotropization. 
Extrapolating our results to conditions expected for 
heavy- ion collisions at the LHC, we found that for the as- 
sumed magnitude of current fluctuations, A = 0.8, com- 
plete isotropization within our framework occurs at ~ 
6.5 fm/c. Further extrapolating our numerical results to 
A = 1.6 which is required in order to achieve occupation 
numbers consistent with quantum fluctuations, one finds 
isotropization times on the order of 5 fm/c. 

We note that although the early indications from ideal 
hydrodynamics would imply that this time scale is much 
too long, in recent years it has emerged that there is 
very little experimental constraint on the degree of local 
momentum-space anisotropy in the quark gluon plasma. 
In HEL the precise time scale for isotropization depends 
on the choice of the amplitude of the initial current fluc- 
tuations and as a consequence the amplitude of the soft 
gauge fields at early times. We have chosen the magni- 
tude of these fluctuations based on studies of the breaking 
of boost invariance in the glasma by quantum fluctua- 
tions. Of course, one can shorten the isotropization time 
by increasing the magnitude of the initial fluctuations 
used; however, within the hard-expanding loop frame- 
work one runs the risk of violating the assumption that 
the energy of the system is dominated by the hard de- 
grees of freedom. Once again this imposes a limit on what 
can be achieved through hard-loop simulations and calls 
for more comprehensive methods to tackle the problem 
which can properly include the back-reaction. 

Our final finding concerned the induced spectrum of 
the unstable soft modes. We found a Boltzmann distri- 
bution of longitudinal energies instead of a power law 



distribution as was found in static simulations. Extrapo- 
lating to RHIC and LHC conditions, this result seems to 
imply that one can achieve longitudinal thermalization of 
the quark gluon plasma on time scales of 1 fm/c. Early 
color glass condensate simulations demonstrated that the 
initial gauge field configurations were transversally ther- 
mal [95] and our results indicate that the system also 
quickly becomes thermal in the longitudinal direction. 

The longitudinal thermalization we see is particular to 
non-Abelian gauge theories. In general, there are two 
effects occurring: (1) mode amplification due to plasma 
instability and (2) mode-mode coupling due to nonlinear 
interactions. In an Abelian plasma only mode amplifica- 
tion occurs and one does not see the emergence of a lon- 
gitudinally thermalized spectrum. One needs the mode- 
mode coupling to spread the deposited energy across 
large ranges of momenta quickly. In the non-Abelian 
case, we have checked different lattice sizes, lattice spac- 
ing, etc. and the rapid emergence of a longitudinally ther- 
malized spectrum seems to be quite robust. 

We note that our finding of an exponential longitudi- 
nal energy spectrum is not in contradiction with the pure 
Yang-Mills simulations of Ref. [82] which found the emer- 
gence of a power-law spectrum, since the power-law spec- 
trum observed therein only emerged at quite late times, 
r > 150 fm/c. At early times Ref. [32] also found what 
appears to be an exponential distribution in the longi- 
tudinal energy spectrum. We also note that usually one 
sees power law spectra energy in the infrared. Due to 
having to use many auxiliary fields we were limited to 
40 2 x 128 lattices. In the future we plan runs on larger 
lattices in order to more carefully determine the infrared 
part of the spectrum. 

Our study, however, is not without caveats. In order 
to have a tractable way to treat the time-dependent hard 
particles, we approximated them as a longitudinally free 
streaming ensemble. This is an extreme assumption that 
should be relaxed, if possible, in the future. Some work 
along this direction has been started in Ref. 99 where 
the authors were able to derive an evolution equation 
for a stable uniform chromoelectric field in an arbitrary 
time-evolving anisotropic background. It would be very 
interesting to see if the method employed in Ref. [95] can 
be extended to the entire stable and unstable mode spec- 
trum. This caveat aside, it's interesting that even with 
such an extreme particle pressure anisotropy being devel- 
oped, the chromo-Weibel instability is able to isotropize 
the system on time scales relevant for heavy-ion colli- 
sions. 

The second important caveat is that we did not in- 
clude the effect of the back reaction of the hard particles 
on the unstable soft gauge fields. Our results seem to 
indicate that the fields grow unabated until there will be 
a significant backreaction. Of course, as soon as the field 
amplitudes become large enough for any back reaction to 
occur, it is possible that this could reduce the anisotropy 
of the hard-particles and reduce the rate of growth of the 
unstable soft modes. 
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In the context of our numerical results, the observation 
of continued unstable mode growth places an upper limit 
on the amount of time over which we can trust our hard- 
loop simulations; however, we find that assuming that 
the initial fraction of the energy carried by soft fields is 
small compared to the hard scale there is a window of 
time over which we can reliably simulate the dynamics. 
Our results indicate a very fast path to isotropization 
within this window of reliability. Addressing the ques- 
tion of the late time dynamics of the system is not pos- 
sible within this frameworkj^] however, our study might 
serve as a benchmark for future simulations that include 
backreaction in an expanding plasma. 

In the future one might use hard-loop simulations to 
study the early time dynamics of the quark gluon plasma 
and the role unstable modes play. One can address in- 
teresting phcnomcnological questions such as measuring 
the shear viscosity due to plasma instabilities and study- 
ing particle transport properties such as energy loss and 
momentum-space diffusion. In addition, the momentum- 
space anisotropy dependence of many important heavy- 
ion collision observables such as jet energy loss, pho- 
ton production, dilepton production, heavy quark energy 
loss, heavy quarkonium suppression etc. have been com- 
puted [61)1 1100f4l 15j . It would, therefore, be interesting 
to study the effect of our time-dependent evolution on 
these observables as a possible signature of the plasma 
instability in heavy ion collisions. 

We note that there are now many groups studying the 
thermalization, isotropization, and anisotropic signatures 
of the quark gluon plasma in the strong coupling limit 
using the AdS/CFT correspondence [29H3T1 ITT6HI28] . It 
would be interesting to compare and contrast the pre- 
dictions for experimental observables coming from the 
weakly-coupled and strongly-coupled frameworks. 

Finally, there have also been some recent studies that 
have suggested that there is an inverse particle number 
cascade leading to Bosc-Einstein condensation of soft- 
gauge fields |129H131j . How a long-lived condensate 
can emerge in a non-Abelian gauge is an open question. 
Based on our results it is hard to judge whether this pos- 
sibility is borne out, since we do not directly obtain the 
particle number spectra but instead the energy spectra. 
Determining the nature of the low momentum number 
spectra is complicated by gauge invariance issues; how- 
ever, measurements of this spectra in fixed-anisotropy 
hard-loop simulations [55J 153 EI] and pure Yang Mills 
with high occupancy [131H133] have so far shown no evi- 
dence of occupation numbers exceeding / ~ l/a s at late 
times. That being said it would be interesting to see if a 



In this paper, we have concentrated on the phenomenology of 
unstable modes in a longitudinally free streaming background. 
For an in-depth analysis of the path to isotropy in the asymp- 
totically small coupling limit, including late time dynamics, we 
refer the reader to Refs. [761177] where parametric estimates have 
been made. 



time-evolving condensate, perhaps in the form of an over- 
populated condensate of plasmons (chromoelectric oscil- 
lations), could play a role in QGP thermalization and 
isotropization. 
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Appendix A: Lattice equations of motion 

In this appendix we introduce the dimensionless lat- 
tice variables we use in simulating the dynamics of the 
soft color fields. We then explicitly write the discretized 
equations of motion and initial conditions used in the 
main body of the paper. In this paper we consider the 
non-Abelian SU(2) group; however, the equations below 
are independent of the gauge group considered. 



1. Lattice variables 

We begin by defining dimensionless lattice variables 
which will be used in the simulation. We introduce three 
lattice spacings: a which is the dimcnsionful transverse 
spatial lattice spacing, e which is the dimcnsionful tem- 
poral lattice spacing, and a v which is the dimensionless 
lattice spacing in the r\ direction. We rescale space and 
time 

x = x/a, y = y/a, f) = V , 

r = r/a, e = e/a, (Al) 

With these definitions we can rescale the field variables, 
conjugate momenta, and currents and introduce lattice 
variables with "hats" 

A 1 = gaA 1 , A v = gA v , 

f[ s ; = gaUi , ft" = ga 2 W , (A2) 

and 

W = aW, J T = a 3 f, 

/ <> '■/ • 3*=a*j*. (A3) 

Finally, we rescale the isotropic Debye mass via rhp = 
amr>. Performing this transformation on the Hamilto- 
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nian density we find 
1 



n 



g 2 a 3 



tr 



F 
j.2 m 



rW + Ft 



(A4) 



In the following subsections we will drop the "hats" on 
symbols. From this point on in this appendix, all vari- 
ables can be assumed to be dimensionless lattice vari- 
ables. 



2. Plaquettes and Staples 

We can translate the continuum equations of motion 
into gauge-invariant lattice equations of motion by using 
standard plaquette and staple operators. For the trans- 
verse coordinate-rapidity plaquettes we have 



iN, 



(F kv ) a = — tr [r a U DM ] 



(A5) 



where k £ {x 7 y} 7 a is a color algebra index, a e 
{l,--- ,iV c 2 -l}, and Ua^ix) = U li (x)U v (x + n)Ut(x + 
v)Ul(x) is a standard lattice plaquette variable with {7 M 
being a parallel transporter in fj, direction 



Ui = cxpi-iA 1 ) , 
U v — exp(?a T? A r) ) . 



(A6) 
(A7) 



Products like F 2 V which appear in the energy density 
can also be expressed in terms of plaquette variables. For 
this application we need two different combinations, F 2 { 



and F 2 y . These are 



tr F l = Jf i 1 ~ Jf c tr t Re u n,*]J > (A8) 



trf£ w = 2(l- — tv[ReU n<xy 



(A9) 



Finally, we can rewrite the necessary covariant deriva- 
tives acting on the field strength tensor as 



(DjjFjjj) 



iN r tr 



iN r 



■ tr 



iN r 



tr 



T a U k (T,x)J2sUr,x) 



T a U n (T, X )J2sl J (T, 



(A10) 



(All) 



(A12) 



where S is the gauge link staple 

SUt, x) = U v {r, x + M )C/t( T; x + v )ul^ x ) . (A13) 



Note that the sums in ( |A12[ ) run over both positive and 
negative directions. 



3. Transformation of the W fields to a compact 
domain 



In order to better describe W in the y (shifted rapidity) 
direction we introduce a velocity-like variable u, — 1 < 
u < 1, defined by 



y = atanh(w) , dy 



1 



l-u 2 



rdu. 



(A14) 



This has the effect of giving more lattice points around 
y = 0, where the W functions are rapidly varying. 

Using sinh 2 (y) = u 2 1 (1 — u 2 ) and cosh 2 (y) 
1 : 1 — u 2 ) we can rewrite (3.331 as 



<9 t W(t, x, n; <i>, u) = -Vl-u 2 d'AW 
-- (D„W - (1 - u 2 )d u W) 



1 



/(r,r iso ,u) 

J2 



-vm, 



=n" + "(i-^)^ 



where 



r iso Vl- a 



t 2 u 2 ^ 2 
/(r, r iso , u) = ( 1 + — 

Tso 1 - u 



(A15) 



(A16) 



The currents are then given by 



,2 flit 



2tt 



9 



2tt 



n 2 /-27T 



2r 



2tt 



j du(l- 


u 2 ) *W(r,x, r);<f> 


tt) , 






(A17) 


j duv l (l 




<M) , 






(A18) 


j duu{\ 


-u 2 )-i W(t,x,77; 


0,u) , 






(A19) 



where, as usual, v l — (cos </>, sin </>) with i € {x, j/}. 

4. Lattice equations of motion 

We will express the equations of motion in terms of 
gauge links U and chromoclcctric fields II. Both t/'s and 
ITs live on links (between sites) so all of their spatial ar- 
guments have an implicit +1/2 shift. In some cases this 
1/2 is made explicit for maximum clarity. Temporally 
ITs also live between sites. The link variables U, how- 
ever, temporally live on sites. The W's and j's live on 
sites both spatially and temporally. We use a lattice with 
N± sites in the x and y directions and N n sites in the r\ 
direction. The fields are assumed to be periodic in all di- 
rections. We use a leapfrog algorithm in which the conju- 
gate momenta are updated first using fixed links/currents 
and then the link variables and W-fields are evolved using 
the updated conjugate momenta [134H138) . 
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The resulting Yang-Mills update equations are 

n »( T + 7>, x ,v) = 11,(1-- ^,x,r?) +re ( j* (r, x, 77) + DjFj^r, x, 7/) + -^£> n F^(r, x, r?) 



IF (r + -, x, 77) = IF (r - x, ry) - - (r 2 ^ vg (r, x, rj) + D z F lv (r, x, rj) 
Ui(T + e,x,r)) = cxp(-ie t _1 IL^t + |,x,r?)J [/,•(•?-, x, 77) , 
t/^(r + e,x,ry) = exp (+ i e t a n IF (r + ^,x,r))j U v (t,x,t)) , 



where 



Javg( r ' X ''7) 
iav S ( r > X ^) 



f(r, x, 77) + [// (r, x, 77) j l (r, x + e;, r))Ui (r, x, 77) 
F(t,x,?7) + ^(r,x,?7)j I '(T,x,?7 + l)C/^(r, x, 77) 



To discretize the W fields we use a rectangular lattice in <fi-u space of size x N u and 

4> n = 2m/ N$ , 
u m = -l + (2m + l)/N u , 

where n g {0, • • • , N4 - 1} and m e {0, • • • , N u - 1}. 



(A20) 

(A21) 
(A22) 
(A23) 

(A24) 
(A25) 

(A26) 
(A27) 



I 

The update equations for the W fields then take the transverse and rapidity directions, respectively 
form 



D 



W(t + e, x, 77; (p, u) — 

W(r -e,x, 77; (f>,u) + 2e\ -^l-u 2 v l D?W 



(D*W - (1 - u 2 )d s u W) 



/(r,r iso ,u) 



-vmr s - - 



"'l-^h'F, 



z?7 



(A28) 



-U v (ri-l)<p{Ti-l)UtiTi-l) 



(A29) 



-Uiixi-lMxi-lplixi-l) 

and is a symmetric derivative in u space 
_ tp(u m+ i) - ip(u m -i) 



(A30) 



d b ip(u m ) 



2Au 



(A31) 



where Au = 2/N u . The averaged conjugate momenta, 



where F iri is computed using plaquettes via Eq. (A5|, n^ vg and II?' appearing in (A28) are averaged both 



Df and D° are symmetric covariant derivatives in the spatially and temporally 



nr g c™)=i e E E E tttu 



(7 X — ±1 (Tu=±l cr„— ±1 cr T — ±1 



n^,x,77)^l E E E E^(r+^,, + ^,77+^,77+^ 



(A32) 
(A33) 



(T x —±1 a v — ±1 cr,,— ±1 cr T — ±1 



where we have indicated explicitly the fact that the El's stands for the parallel transporter necessary to bring the 
live on links (halfway between sites) for clarity and Vj- 
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conjugate momenta to the same site. 

The currents are computed from the W fields via 



D 



N^N U 



" L D 



(A34) 

^'(l-u^Wfox^.u), 

n,m 

(A35) 



f(r, x, r,) = -p^- £ u (1 - u 2 )- § W(r, x, r?; 0, «) , 

^ u n,m 

(A36) 

and we monitor Gauss' Law by periodically checking 



tr 



kv- Yl ( T ' X ' + D V n avg (T, X, 7?) 



iV 2 A^ 



£fnr*(T,x,»j) 



(A37) 



We compute the discretized transverse and longitudinal 
contributions to the field energy density £ via 

^ = ^£tr[T- 2 F 2 t+ r- 2 n 2 ], (A38) 



?2 r,r,A f t- P2 



(A39) 



X. 7} 

where tri* 1 ^ and tiF^ y are computed using Eqs. (A8) 



and (A9) 



Appendix B: Choice of lattice parameters 

In this appendix we detail the constraints which should 
be obeyed in order for our simulations to properly de- 
scribe the soft gauge field dynamics. Since the soft scale 
is time dependent, we have to choose parameters which 
allow for a faithful representation of the infrared and ul- 
traviolet physics during the entirety of the simulation. 

The physical (dimensionful) parameter m 2 D is the De- 
bye mass at time Ti so . In terms of the gluon liberation 
factor c which is 0(1) (— 1.1 according to Lappi [9"4"] . 
c = 2 In 2 w 1.386 according to Kovchegov [91]) one has 



m% r iso r « 0.93c(<3 s t ) . 



(Bl) 



In the text we use Q s tq = 1 and c = 2 In 2 from 
Kovchegov [51] . This gives m 2 D T lso TQ = 1.285. For large 
anisotropy one finds 



2 I \ 71 2 / 



which can be taken as the typical (time-dependent) soft 
momentum scale. With our choice of c = 2 In 2, we 



have m^r) ss 1.0 [tqt) 1 / 2 . For RHIC energies one has 
t q 1 = Qs ~ 1-4 GeV and at current LHC energies one 
has Q s ~ 2 GeV. At RHIC and LHC energies t = Q7 1 
corresponds to 0.14 fm/c and 0.1 fm/c, respectively. 

On a lattice with periodic boundary conditions, the 
size of the lattice determines the infrared cutoff in full 
wavelengths and the lattice spacing determines the ul- 
traviolet cutoff via the smallest half-wavelength. In our 
expanding system, the transverse UV cutoff is constant 
in time and given by n/a, whereas the soft momentum 
scale is decreasing in time. It is therefore sufficient to 
ensure 



'max = - > TOoo(To) 



(B3) 



so that we should demand a < Itq. 

In longitudinal direction, the effective UV cutoff is de- 
creasing in time according to n/(Ta v ). We may choose 
for example 

— — a n ~ a, (B4) 



to have comparable transverse and longitudinal UV cut- 
offs in an average sense. More importantly, the maximal 
longitudinal wave number of unstable modes increases in 
time, so ^ max should be a large number, 



= — > 30 . 



(B5) 



Since the hard-expanding-loop framework is designed 
to treat the soft sector of the dynamics, it is somewhat 
more important to properly treat the infrared scale. In 
the longitudinal direction v m i n = 2i: / {N^a^) should be 
made as low as possible. There are no important unstable 
modes with v much smaller than 5, but v m in also sets the 
spacing between mode numbers. We should therefore aim 
at 



2tt 



<5. 



(B6) 



In the transverse direction, the semianalytic results [HH 
l86| suggest that we should have 

fcmin-v^«0.2r - 1 . (B7) 



As our canonical set of parameters in the results sec- 
tion we use N T = 40, N v = 128, a v = 0.025, a = Qj 1 , 
and To = Qj 1 ■ Checking the transverse infrared cutoff 
one finds /c m i n = 0.157 Q s < 0.2 Q s as required. Check- 
ing the transverse ultraviolet cutoff one finds fc ma x = 
(B2) ttQs > 1.005 Q s as required. Checking the longitudi- 
nal infrared cutoff one finds f m in = 1-96 < 5 as required. 
Finally, checking the longitudinal ultraviolet cutoff one 



finds v„ 



125.7 > 30. 
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Case 



Std. Dev. (o) 



Transverse Vector Potential (Ai 



gA/a 



Longitudinal Vector Potential (A v ) gA/(aa 



1/2 
) 

1/2, 



Transverse Conjugate Momentum (IT) gA/a^ 2 
Longitudinal Conjugate Momentum (IT,,) gaA/a^ 2 
Current fluctuations (W) 



A/ a y 2 



TABLE I. Transverse and longitudinal lattice spacing scaling 
for a variety of different initial condition types. 



Appendix C: Initial Conditions 

In this appendix we collect details of the initial condi- 
tions used in the simulations and some information about 
lattice initial conditions in general. 



1. Gaussian random variables 

We now discuss the scalings necessary when sampling 
lattice variables from Gaussian distributions. For com- 
pleteness we list all possible types of initial conditions; 
however, in the body of the text we use exclusively ini- 
tial conditions based on current fluctuations. We then 
give some more details about the precise implementation 
of the current fluctuation initial conditions used in the 
body of the text. 

It is common to use uncorrelated Gaussian random 
noise as the initial condition for either fields or current 
fluctuations. In the case of uncorrelated transverse vector 
potentials, for example, one assumes that in the contin- 
uum limit 

(A?(T ,x 1 ,rn)A b j (T ,X2,m)) = 

A 2 S ab S^ 2 \ Xl ~ x 2 )%! - rft) , (CI) 

where x^ = (x, y) is a purely transverse two-vector. In 
order to translate this statement into something useful 
for the lattice initial conditions we should convert to di- 
mensionless variables on the left and right hand sides. 
In doing so we make use of the rescalings specified in 
Eqs. (All and (A2) and the Dirac delta function iden- 
tity S (ax) — <5(x)/|a| to obtain (in terms of the lattice 
variables introduced in Eqs. (A2) and (A3)) 



(C2) 



In practice, this means that the Ai variables should be 
Gaussian random numbers with a standard deviation of 
(j = gAja\j 2 . 

Using similar arguments we can derive the following 
lattice correlation functions in the case that we initialize 



5__^_raf)r £ 



q 2 A 2 



longitudinal vector potentials 

or transverse momenta 
(n^To.x^^n^ro.x^a)) = 

or longitudinal momenta 

(n^x^on^x^)) = 

or auxiliary fields 



— 6 6 ^ x ^X2^Viy2^VlV2^<t>l<t>2 ■ 



(C3) 



5 SijS x ± x x8 r i iri2 , 
(C4) 



n 2 n 2 A 2 

9 a A S ab S S 

Cinq 



(C5) 



(C6) 



To summarize, when using Gaussian random initial con- 
ditions on anisotropic lattices, one should choose the 
standard deviations shown in Table |T| Moreover, unless 
initial fluctuations are only set up for the gauge fields Af 
and A i, a projection to satisfy the Gauss law constraint 
(3.26) is needed. In our simulations we have however 



used a different setup which we now discuss. 



2. Initial Condition Setup 

The analytic study of collective modes in anisotrop- 
ically expanding ultarelativistic plasmas |86j has found 
that the initial fluctuations in (only) induced currents 
versus only initial fluctuations in collective fields reduces 
considerably the delay of the onset of the plasma instabil- 



ities. As discussed in Sect. IV this means that such ini 



tial conditions dominate over all other possibilities, and 
it is therefore sufficient to concentrate on initial fluctua- 
tions in the W fields which directly encode the induced 
currents. 



a. Longitudinal Current Initial Conditions 

For oblate anisotropy, fluctuations in longitudinal cur- 
rents give rise to stable plasmon modes, and in the 
Abelian case they do not lead to any plasma instabili- 
ties. We have used this to test our code for unphysical 
instabilities (see App. [d|. 

The simplest initial fluctuations consistent with Gauss' 
law which achieve this are fluctuations in only the W ar) 
components that are independent of <fi an d y (thereby 
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FIG. 12. (Color online) Collected numerical tests of unstable mode growth. Each subpanel shows the chromofield total energy 
density evolution subject to variation of various parameters. The subcaptions contain a description of the parameters which 
are varied. 
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ensuring that initially j T — 0) but nothing else 

(IF a "(r , m; 0i, j/i)W 6 "(r , x£, 772; 2 , y 2 )) 

A 2 



IF" = 0, £7. , 1 = 1 



"aj^a;^- ^!^) 



(C7) 



Transversal Current Initial Conditions 



In order to provide seed fields for Weibel instabilities, 
longitudinal current fluctuations do not play an impor- 
tant role (for oblate anisotropies). For simplicity we have 
therefore only considered transverse current fluctuations 
by only initializing W a% fields. Because we have used 
rather fine lattices in the r\ direction, Gaussian random 
noise would correspond to very high UV noise in longi- 
tudinal wave numbers even beyond the scale which sepa- 
rates soft and hard modes, while hard modes are already 
integrated out. We have therefore introduced a mode 
number cutoff A„ such that v max — A„i/ m i n with v m [ n — 
2tt / (Nrja^) . Again, the simplest initial fluctuations con- 
sistent with the Gauss law are obtained by requiring that 
the W al components are independent of tj> and y, and 
thus initially j T = 0, while setting all other fields to zero 
initially. This is now done in terms of the Fourier compo- 
nents W al with W al {. . . , 77, . . .) = Y,u W ai (- ■ ■ , v, . . .)e iur > 
according to 

(W al (to , xt, vi ; 0i , yi ) W bj (r , x^ , u 2 ; <j> 2 , y 2 ) ) = 
A 2 h 



6 a 5 v S x x x xS Ul - U2 9(u max - \ui\), 



0, u s+i = i Nc 



- - 0. 



(C8) 



Appendix D: Numerical Tests 



In this section we collect various numerical tests such 
as varying the lattice spacing, lattice size, spectral cut- 
offs, and velocity-space resolution. In Fig. [12] we collect 
six different tests. The variation with the random seed 
used for generating the necessary pseudorandom num- 



bers used in the initial conditions is shown in Fig 12(a) 
As we can see from this figure there is a fair amount of 
variation with the random seed used; however, the results 
are all qualitatively the same. In the results section our 
main results are averaged over the set of runs shown in 



Fig 12(a) 



The variation with the ultraviolet longitudinal mode 
cutoff used for initializing the initial current fluctuations 
via the auxiliary VV fields is shown in Fig |12(bj| As can 
be seen from this figure there is a rapid convergence as 
the ultraviolet cutoff, A. v i' m i n , is increased. The set of 
runs shown in the main body of the text uses A, y = 8. 



The variation with the transverse lattice spacing while 
holding the transverse lattice size fixed is shown in 
Fig 12(c) This represents a test of the approach to 



the continuum as the transverse lattice resolution is in- 
creased. In the transverse plane we sample Gaussian 
random numbers which means as the lattice spacing de- 
creases the transverse configurations will be dominated 
by the high transverse momentum part of the fluctua- 
tions. This is evidenced by the fact that the initial en- 
ergy density deposited in the fields by the current fluctu- 
ations increases rapidly as one approaches the transverse 
continuum limit. One could remove this artifact by im- 
plementing a transverse mode cutoff on the lattice, but 
at this point in time we have not yet done so. In the 
results section our standard set of runs uses Nt = 40. 

The variation with the transverse lattice size while 
holding the transverse lattice spacing fixed is shown in 
Fig 12(d) In this case we see a rather large effect. In the 
limit that Nt — > 1 while holding a fixed, one approaches 
a one-dimensional system which exhibits a faster growth 
rate due to less mode competition. We have verified that 
in this limit we reproduce our previously obtained re- 
sults from Ref. [55]. The faster growth seen compared 
to Ref. [55] is due to the use of the more general initial 
conditions which include current fluctuations [86j . In the 
results section our standard set of runs uses Nt = 40. 
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FIG. 13. (Color online) Evolution of a stable configuration ini- 
tialized with Abelian longitudinal currents for different sized 
velocity lattices. 

The variation with the longitudinal lattice spacing 
while holding the longitudinal lattice size fixed is shown 
in Fig 12(e) Due to the fact that we have implemented 



an ultraviolet cutoff on fluctuations in the ?7-direction, 
we see a very nice convergence as the lattice resolution 
in the r\ direction is increased. In the results section our 
standard set of runs uses Nt = 40. In the results section 
our standard set of runs uses N v = 128. 

The variation with the longitudinal lattice size while 
holding the longitudinal lattice spacing fixed is shown in 
Fig |12(f)j Once again we see only small variation with 
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the assumed longitudinal lattice size, with the late-time 
variations being consistent with those coming from ran- 
dom seed variation. In the results section our standard 
set of runs uses N„ = 128. 



Finally, in Fig. 13 we show the evolution of a stable 
Abelian configuration initialized with Abelian longitudi- 
nal currents for various different velocity lattice resolu- 
tions N u x e {64 x 16, 128 x 16, 128 x 32, 128 x 48}. 
For this simulation the lattice spatial size was x N v = 
32 2 x 32 with transverse lattice spacing of a = 0.1 fm and 
longitudinal lattice spacing of a„ =0.025. The initial time 
was taken to be tq = 0.1 fm/c and we used T[ SO = 0.01 



fm/c. For the temporal time step we use e = 10 -3 fm/c. 
With these initial conditions, the field energy should de- 
cay steadily after the initial peak. This test turns out to 
be very sensitive to the velocity-space resolution, i.e. the 
number of W fields. If this resolution is too crude, the 
field energy even grows at late times. Fig. [13] shows that 
with a velocity lattice size of N u x N$ = 128 x 32 there is 
already good convergence to the correct time evolution 
of the system. Unstable modes are in fact less sensitive 
to the velocity resolution in the <fi direction; however, 
being cautious we have performed all simulations using 
N u x = 128 x 32. 
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